
data_1 <-
  data_1 %>%
  mutate(q15_java_binary = case_when(q15_survey_test_java %in% c(3, 4) ~ 1,
                                     q15_survey_test_java %in% c(1, 2) ~ 0,
                                     TRUE ~ NA_real_),
         q15_muslim_binary = case_when(q15_survey_test_muslim %in% c(3, 4) ~ 1,
                                       q15_survey_test_muslim %in% c(1, 2) ~ 0,
                                       TRUE ~ NA_real_)
  )


############################################
#TESTING EFFECT#############################
############################################

test_mod1a <- lm(java_index_skd_bw_1 ~ fail_skd*q15_java_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod1b <- lm(java_index_skd_bw_1 ~ fail_skd*q15_muslim_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod2a <- lm(non_java_index_skd_bw_1 ~ fail_skd*q15_java_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod2b <- lm(non_java_index_skd_bw_1 ~ fail_skd*q15_muslim_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod3a <- lm(regional_index_skd_bw_1 ~ fail_skd*q15_java_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod3b <- lm(regional_index_skd_bw_1 ~ fail_skd*q15_muslim_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod4a <- lm(religious_index_skd_bw_1 ~ fail_skd*q15_java_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod4b <- lm(religious_index_skd_bw_1 ~ fail_skd*q15_muslim_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod5a <- lm(corruption_index_skd_bw_1 ~ fail_skd*q15_java_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod5b <- lm(corruption_index_skd_bw_1 ~ fail_skd*q15_muslim_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))

test_mod6a <- lm(national_index_skd_bw_1 ~ fail_skd*q15_java_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))
test_mod6b <- lm(national_index_skd_bw_1 ~ fail_skd*q15_muslim_binary, data = data_1 %>% filter(abs(forcing_skd) < 5))

observations <- c(nobs(test_mod1a), nobs(test_mod2a), 
                  nobs(test_mod3a), nobs(test_mod4a), 
                  nobs(test_mod5a), nobs(test_mod6a))

test_mod1a <- coeftest(test_mod1a, vcov=vcovHC(test_mod1a,type="HC0"))
test_mod1b <- coeftest(test_mod1b, vcov=vcovHC(test_mod1b,type="HC0"))

test_mod2a <- coeftest(test_mod2a, vcov=vcovHC(test_mod2a,type="HC0"))
test_mod2b <- coeftest(test_mod2b, vcov=vcovHC(test_mod2b,type="HC0"))

test_mod3a <- coeftest(test_mod3a, vcov=vcovHC(test_mod3a,type="HC0"))
test_mod3b <- coeftest(test_mod3b, vcov=vcovHC(test_mod3b,type="HC0"))

test_mod4a <- coeftest(test_mod4a, vcov=vcovHC(test_mod4a,type="HC0"))
test_mod4b <- coeftest(test_mod4b, vcov=vcovHC(test_mod4b,type="HC0"))

test_mod5a <- coeftest(test_mod5a, vcov=vcovHC(test_mod5a,type="HC0"))
test_mod5b <- coeftest(test_mod5b, vcov=vcovHC(test_mod5b,type="HC0"))

test_mod6a <- coeftest(test_mod6a, vcov=vcovHC(test_mod6a,type="HC0"))
test_mod6b <- coeftest(test_mod6b, vcov=vcovHC(test_mod6b,type="HC0"))


#Javanese advantage table

table <- list(test_mod1a, test_mod2a, test_mod3a, test_mod4a, test_mod5a, test_mod6a)

note_text <- paste(" Beta coefficients from OLS regression with interaction terms. Standard errors were calculated using the Huber-White (HC0) correction. 
                   The outcomes measure are indexed values capturing (1) Javanese preferentialism among Javans and (2) among non-Javans, (3) regional preferentialism,
                   (4) religious resentment, (5) perceptions of corruption, (6) national identification.")

table = stargazer(table, 
                  type = 'latex', 
                  title = "The Effect of Basic Competence Examination (SKD) Failure, By Perceptions of Javan Advantage",
                  label = 'tab:testing_effect_javan_hte',
                  model.names = F,
                  model.numbers = T,
                  digits = 2,
                  column.separate = c(2, 1, 1, 1, 1),
                  column.labels = c("Java. Pref.", "Reg. Pref.", "Relg. Resent.", "Corrup. Percep.", "Natl. ID"),
                  multicolumn = T,
                  dep.var.labels = NULL, 
                  add.lines = list(c("Subset", "Javan", "non-Javan", "---", "---", "---", "---"),
                                   c("Observations", observations),
                                   c('Bandwidth', rep(c('1\\%'), 18))),
                  covariate.labels = c("Failed SKD", "Believe Javans Advantaged", "Believe Javans Advantaged X Failed SKD"),
                  #star.cutoffs = c(0.000000000000000000000000000000000000000000000000000000000000),
                  #float.env = 'sidewaystable',
                  keep.stat = c("n"),
                  notes = NULL,
                  notes.align = 'l')


write_latex(table[-c(10, 11, 12, 18, 21, 24, 27, 34)], note_text, './_4_outputs/tables/table_a11.tex', .8)



#muslim advantage table

table <- list(test_mod1b, test_mod2b, test_mod3b, test_mod4b, test_mod5b, test_mod6b)


observations <- c(nobs(test_mod1b), nobs(test_mod2b), 
                  nobs(test_mod3b), nobs(test_mod4b), 
                  nobs(test_mod5b), nobs(test_mod6b))

note_text <- paste(" Beta coefficients from OLS regression with interaction terms. Standard errors were calculated using the Huber-White (HC0) correction. 
                   The outcomes measure are indexed values capturing (1) Javanese preferentialism among Javans and (2) among non-Javans, (3) regional preferentialism,
                   (4) religious resentment, (5) perceptions of corruption, (6) national identification.")

table = stargazer(table, 
                  type = 'latex', 
                  title = "The Effect of Basic Competence Examination (SKD) Failure, By Perceptions of Muslim Advantage",
                  label = 'tab:testing_effect_muslim_hte',
                  model.names = F,
                  model.numbers = T,
                  digits = 2,
                  column.separate = c(2, 1, 1, 1, 1),
                  column.labels = c("Java. Pref.", "Reg. Pref.", "Relg. Resent.", "Corrup. Percep.", "Natl. ID"),
                  multicolumn = T,
                  dep.var.labels = NULL, 
                  add.lines = list(c("Subset", "Javan", "non-Javan", "---", "---", "---", "---"),
                                   c("Observations", observations),
                                   c('Bandwidth', rep(c('1\\%'), 18))),
                  covariate.labels = c("Failed SKD", "Believe Muslims Advantaged", "Believe Muslims Advantaged X Failed SKD"),
                  #star.cutoffs = c(0.000000000000000000000000000000000000000000000000000000000000),
                  #float.env = 'sidewaystable',
                  keep.stat = c("n"),
                  notes = NULL,
                  notes.align = 'l')


write_latex(table[-c(10, 11, 12, 18, 21, 24, 27, 34)], note_text, './_4_outputs/tables/table_a12.tex', .8)

